Chapter 1: Predictive modeling

Get prepared

install.packages(c("tableone", "DALEX", "ggplot2", "partykit", "mlr3", "mlr3learners", "ranger", "mlr3tuning", "paradox"))

Part 1: Introduction to predictive modeling + EDA

The purpose of this tutorial is to present a life cycle of a single predictive model for a problem related to binary classification (but I deliberately don’t mention logistic regression). Along the way, we will tackle various interesting topics such as model training, model verification, hyperparameter tuning, and exploratory model analysis.

The examples presented here are inspired by three textbooks books: The Elements of Statistical Learning with mathematical foundations, The mlr3 book, presenting software package mlr3 designed for advanced model programming, Explanatory Model Analysis overviews methods for model exploration and visualisation. Note that responsible modelling requires knowledge that cannot be learned in two days.
So, after this introduction, I highly recommend checking out these books.

Why should I care?

Predictive models have been used throughout entire human history. Priests in Egypt were predicting when the flood of the Nile or a solar eclipse would come. Developments in statistics, increasing the availability of datasets, and increasing computing power allow predictive models to be built faster and faster.

Today, predictive models are used virtually everywhere. Planning the supply chain for a large corporation, recommending lunch or a movie for the evening, or predicting traffic jams in a city. Newspapers are full of interesting applications.

But how are such predictive models developed? In the following sections, we will go through a life cycle of a predictive model. From the concept phase, through design, training, checking, to the deployment. For this example, we will use the data set on the risk of death for Covid-19 patients after SARS-COV-2 infection. But keep in mind that the data presented here is artificial. It is generated to mirror relations in real data, but do not contain real observations for real patients. Still, it should be an interesting use-case to discuss a typical lifetime of a predictive model.

Tools

These materials are based on two R packages: mlr3 for model training and DALEX for model visualization and explanation. But there are more packages with similar functionalities, for modelling other popular choices are mlr, tidymodels and caret while for the model explanation you will find lots of interesting features in flashlight and iml.

The problem

The life cycle of a predictive model begins with a well-defined problem. In this example, we are looking for a model that assesses the risk of death after diagnosed covid. We don’t want to guess who will survive and who won’t. We want to construct a score that allows us to sort patients by risk of death.

Why do we need such a model? It could have many applications! Those at higher risk of death could be given more protection, such as providing them with pulse oximeters or preferentially vaccinating them.

Load packages

library("tableone")
library("DALEX")
library("ggplot2")
library("partykit")
library("mlr3")
library("mlr3learners")
library("ranger")
library("mlr3tuning")
library("paradox")

set.seed(1313)

Conception

Before we build any model, even before we touch any data we should first determine for what purpose we will build a predictive model.

It is very important to define the objective before we sit down to programming, because later it is easy to get lost in setting function parameters and dealing with all these details that we need to do. It is easy to lose sight of the long-term goal.

So, first: Define the objective.

For the purpose of these exercises, I have selected data on the covid pandemic. Imagine that we want to determine the order of vaccination. In this example we want to create a predictive model that assesses individual risks because we would like to rank patients according to their risk.

To get a model that gives a best ranking we will use the AUC measure to evaluate model performance. What exactly the AUC is I’ll talk about a little later, right now the key thing is that we’re interested in ranking of patients based on their risk score.

Read the data

To build a model we need good data. In Machine Learning, the word good means a large amount of representative data. Collecting representative data is not easy and often requires designing an appropriate experiment.

The best possible scenario is that one can design and run an experiment to collect the necessary data. In less comfortable situations, we look for “natural experiments,” i.e., data that have been collected for another purpose but that can be used to build a model. Here we will use the data= collected through epidemiological interviews. There will be a lot of data points and it should be fairly representative, although unfortunately it only involves symptomatic patients who are tested positive for SARS-COV-2.

For the purposes of this exercise, I have prepared two sets of characteristics of patients infected with covid. It is important to note that these are not real patient data. This is simulated data, generated to have relationships consistent with real data (obtained from NIH), but the data itself is not real. Fortunately, they are sufficient for the purposes of our exercise.

The data is divided into two sets covid_spring and covid_summer. The first is acquired in spring 2020 and will be used as training data while the second dataset is acquired in summer and will be used for validation. In machine learning, model validation is performed on a separate data set. This controls the risk of overfitting an elastic model to the data. If we do not have a separate set then it is generated using cross-validation, out of sample or out of time techniques.

covid_spring.csv corresponds to covid mortality data from spring 2020. We will use this data for model training.
covid_summer.csv corresponds to covid mortality data from summer 2020. We will use this data for model validation.

covid_spring <- read.table("data/covid_spring.csv", sep =";", header = TRUE, stringsAsFactors = TRUE)
covid_summer <- read.table("data/covid_summer.csv", sep =";", header = TRUE, stringsAsFactors = TRUE)

Explore the data

Before we start any serious modeling, it is worth looking at the data first. To do this, we will do a simple EDA. In R there are many tools to do data exploration, I value packages that support so called table one.

library("tableone")

table1 <- CreateTableOne(vars = colnames(covid_spring)[1:11],
                         data = covid_spring,
                         strata = "Death")
print(table1)

##                                    Stratified by Death
##                                     No            Yes           p      test
##   n                                  9487           513                    
##   Gender = Male (%)                  4554 (48.0)    271 (52.8)   0.037     
##   Age (mean (SD))                   44.19 (18.32) 74.44 (13.27) <0.001     
##   Cardiovascular.Diseases = Yes (%)   839 ( 8.8)    273 (53.2)  <0.001     
##   Diabetes = Yes (%)                  260 ( 2.7)     78 (15.2)  <0.001     
##   Neurological.Diseases = Yes (%)     127 ( 1.3)     57 (11.1)  <0.001     
##   Kidney.Diseases = Yes (%)           111 ( 1.2)     62 (12.1)  <0.001     
##   Cancer = Yes (%)                    158 ( 1.7)     68 (13.3)  <0.001     
##   Hospitalization = Yes (%)          2344 (24.7)    481 (93.8)  <0.001     
##   Fever = Yes (%)                    3314 (34.9)    335 (65.3)  <0.001     
##   Cough = Yes (%)                    3062 (32.3)    253 (49.3)  <0.001     
##   Weakness = Yes (%)                 2282 (24.1)    196 (38.2)  <0.001

During modeling, exploration often takes the most time. In this case, we will limit ourselves to some simple graphs.

ggplot(covid_spring, aes(Age)) +
  geom_histogram() +
  ggtitle("Histogram of age")

ggplot(covid_spring, aes(Age, fill = Death)) +
  geom_histogram() +
  ggtitle("Histogram of age")

ggplot(covid_spring, aes(Age, fill = Death)) +
  geom_histogram(color = "white") +
  ggtitle("Histogram of age") + 
  DALEX::theme_ema() +
  scale_fill_manual("", values = c("grey", "red3"))

library(ggmosaic)
ggplot(data = covid_spring) + 
  geom_mosaic(aes(x=product(Diabetes), fill = Death)) + 
  DALEX::theme_ema() +
  scale_fill_manual("", values = c("grey", "red3"))

library("pheatmap")
pheatmap((covid_spring[,3:11] == "Yes") + 0)

Transform the data

One of the most important rules to remember when building a predictive model is: Do not condition on future!

Variables like Hospitalization or Cough are not good predictors, beacuse they are not known in advance.

covid_spring <- covid_spring[,c("Gender", "Age", "Cardiovascular.Diseases", "Diabetes",
               "Neurological.Diseases", "Kidney.Diseases", "Cancer",
               "Death")]
covid_summer <- covid_summer[,c("Gender", "Age", "Cardiovascular.Diseases", "Diabetes",
               "Neurological.Diseases", "Kidney.Diseases", "Cancer",
               "Death")]

Your turn

Plot Age distribution for covid_spring and covid_summer.
Calculate tableone for covid_spring and covid_summer.
(extra) In the DALEX package you will find titanic_imputed dataset. Calculate tableone for this dataset.

Part 2: Hello model! First predictive model + How to measure performance

We will think of a predictive model as a function that computes a certain prediction for certain input data. Usually, such a function is built automatically based on the data. But technically the model can be any function defined in any way. The first model will be based on statistics collected by the CDC (CDC stands for Centers for Disease Control and Prevention. You will find a set of useful statistics related to Covid mortality on this page) that determine mortality in different age groups.

In many cases, you do not need data to create a model. Just google some information about the problem.

It turns out that CDC has some decent statistics about age-related mortality. These statistics will suffice as a first approximation of our model.

https://www.cdc.gov/coronavirus/2019-ncov/covid-data/investigations-discovery/hospitalization-death-by-age.html

Lesson 1: Often you don’t need individual data to build a good model.

Create a model

What is a predictive model? We will think of it as a function that takes a set of numbers as input and returns a single number as the result - the score.

cdc_risk <- function(x, base_risk = 0.00003) {
  multip <- rep(7900, nrow(x))
  multip[which(x$Age < 84.5)] <- 2800
  multip[which(x$Age < 74.5)] <- 1100
  multip[which(x$Age < 64.5)] <- 400
  multip[which(x$Age < 49.5)] <- 130
  multip[which(x$Age < 39.5)] <- 45
  multip[which(x$Age < 29.5)] <- 15
  multip[which(x$Age < 17.5)] <- 1
  multip[which(x$Age < 4.5)]  <- 2
  multip * base_risk
}
x <- data.frame(Age = 25, Hospitalisation = "Yes")
cdc_risk(x)

## [1] 0.00045

library("DALEX")
model_cdc <-  DALEX::explain(cdc_risk,
                   predict_function = function(m, x) m(x),
                   type  = "classification",
                   label = "CDC")

## Preparation of a new explainer is initiated
##   -> model label       :  CDC 
##   -> no data available! (  WARNING  )
##   -> target variable   :  not specified! (  WARNING  )
##   -> predict function  :  function(m, x) m(x) 
##   -> predicted values  :  No value for predict function target column. (  default  )
##   -> model_info        :  package Model of class: function package unrecognized , ver. Unknown , task regression (  default  ) 
##   -> model_info        :  type set to  classification 
##   -> model_info        :  Model info detected classification task but 'y' is a NULL .  (  WARNING  )
##   -> model_info        :  By deafult classification tasks supports only numercical 'y' parameter. 
##   -> model_info        :  Consider changing to numerical vector with 0 and 1 values.
##   -> model_info        :  Otherwise I will not be able to calculate residuals or loss function.
##   -> residual function :  difference between y and yhat (  default  )
##   A new explainer has been created!

predict(model_cdc, x)

## [1] 0.00045

The same function can be written in a slightly more compact form as (now it works on many rows)

cdc_risk <- function(x, base_risk = 0.00003) {
  bin <- cut(x$Age, c(-Inf, 4.5, 17.5, 29.5, 39.5, 49.5, 64.5, 74.5, 84.5, Inf))
  relative_risk <- c(2, 1, 15, 45, 130, 400, 1100, 2800, 7900)[as.numeric(bin)] 
  relative_risk * base_risk
}

# check it
x <- data.frame(Age = c(25,45,85))
cdc_risk(x)

## [1] 0.00045 0.00390 0.23700

summary(cdc_risk(covid_spring))

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00003 0.00135 0.00390 0.01849 0.01200 0.23700

table(Death = covid_spring$Death, 
      Prediction.above.005 = cdc_risk(covid_spring) > 0.05)

##      Prediction.above.005
## Death FALSE TRUE
##   No   8946  541
##   Yes   237  276

Wrap the model

In R, we have many tools for creating models. The problem with them is that these tools are created by different people and return results in different structures. So in order to work uniformly with the models we need to package the model in such a way that it has a uniform interface.

Different models have different APIs.

But you need One API to Rule Them All!

The DALEX library provides a unified architecture to explore and validate models using different analytical methods.

Introduction to Responsible Machine Learning

With mlr3 and DALEX

Przemyslaw Biecek and Anna Kozak

07.07.2021

Chapter 1: Predictive modeling

Part 1: Introduction to predictive modeling + EDA

Load packages

Conception

Read the data

Explore the data

Transform the data

Your turn

Part 2: Hello model! First predictive model + How to measure performance

Create a model

Wrap the model

Model performance

ROC

LIFT

Your turn

Part 3: Basics of decision tree and random forest

Test your model

Your turn

Plant a forest

Test your model

Your turn

Part 4: Hyperparameter optimization + Wrap-up

Automated Hyperparameter Optimisation

Define the search space

Set-up the tuner

Tune

Test your model

Sum up

Your turn

Chapter 2: Model exploration

Part 1: Model level analysis - variable importance

Your turn

Part 2: Model level analysis - variable profile

Your turn

Part 3: Instance level analysis - variable attributions

Your turn

Part 4: Instance level analysis - variable profile + Wrap-up

Your turn

Extras

Session info